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Abstract 

This paper concerns with nonuniform sampling and interpolation methods combined with 
variational models for the solution of a generalized image inpainting problem and the restoration 
of digital signals. In particular, we discuss the problem of reconstructing a digital signal/image 
from very few, sparse, and complete information and a substantial incomplete information, which 
will be assumed as the result of a nonlinear distortion. As a typical and inspiring example, we 
illustrate the concrete problem of the color restoration of a destroyed art fresco from its few 
known fragments and some gray picture taken prior to the damage. Numerical implementations 
are included together with several examples and numerical results to illustrate the proposed 
method. The numerical experience suggests furthermore that a particular system of coupled 
Hamilton- Jacobi equations is well-posed. 

AMS subject classification: 35A15, 65M06, 65M32, 68U10, 70H20, 94A08, 94A14, 94A20 

Key Words: signal and image dynamic processing, inpainting, art restoration, variational calculus, 
nonuniform sampling. 

1 Introduction 

Imagine an important and huge art fresco, maybe dated to 1450, right at the very beginning of the 
Italian Renaissance, and imagine that one day a dramatic and catastrophic event destroyed it into 
fragments, maybe during the Second World War. It would be a very hard problem to puzzle up this 
fragmented opera, especially if the original dimensions were large and the number of fragments was 
of several thousands. 

On 11 th March 1944, a group of bombs launched from an Allied airplane hit the famous Italian 
Eremitani's Church in Padua, destroying it together with the inestimable frescoes of Andrea Man- 
tegna et al. contained in the Ovetari Chapel. People picked up and collected what was remaining 
of one of the most important operas of the Italian Renaissance, saving the fragments in some im- 
provised boxes constructed with the wood of the sits of the Church. Several attempts to restore 
these fragments by traditional methods have been done in the last 60 years, without much success. 
Details on "the state of the art" can be found in the booklet 1201 ■ Since 1998, the author has been 
involved in the fascinating attempt to recall to life these frescoes, by using mathematical methods 
and computer based techniques. Recently a fast, robust, and efficient pattern recognition algorithm 
has been developed |19l 121) in order to detect the right position and orientation of the fragments, 
by means of the comparison with an old gray level image of the fresco prior to the damage 1 . This 
method showed to be very effective and, for example, Figure 1 illustrates some detected fragments 
positioned on their original place. Unfortunately the surface covered by the original fragments is 

1 The gray images of the fresco are dated to 1920. 
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Figure 1: Fragmented A. Mantegna's frescoes (1452) by a bombing in the Second World War. 
Computer based reconstruction by using efficient pattern matching techniques 

only 77 m 2 , while the original area was of several hundreds. This means that what we can currently 
reconstruct is just a fraction of what was this inestimable artwork. In particular, it is not known 
what can be the original color of those missing parts. There exist some color images of parts of 
the frescoes, but it is already proved that such colors are not faithful at all and the quality of such 
pictures is very low. So, natural questions raise: Is it possible to estimate mathematically what 
can be the original colors of the frescoes by using the known fragments' information? And, how 
faithful is this estimation? In this paper we want to discuss these questions and to illustrate some 
possible solutions. Clearly, this problem can be reformulated in very different ways for different 
other situations. Therefore, the ideas developed in this context can be put into more general frames 
as we will illustrate in the following. 

The reconstruction of a signal from sparse and nonuniform sampling data is a well known problem 
in information theory [HI d ■ Mathematical methods and numerical algorithms have been developed 
to compute missing parts of signals from few and sparse known sampling information and we refer, 
for example, to classical works of Feichtinger et al. ^] El QUI an d Aldroubi/Grochenig for 
major details. These methods are essentially based on a (quasi-) interpolation of the signal/image by 
means of suitable series expansions of irregularly shifted (translated) basic functions or, in its discrete 
version [1511161 IS"), by series expansion of complex exponentials. Therefore, we may say that relevant 
mathematical methods and concepts maybe useful in image restoration problems are interpolation 
and Fourier /numerical harmonic analysis techniques. Image interpolation can be called, with a 
more "artistic" synonym, "inpainting" . In the last years, image processing problems attracted the 
attention of the mathematical community, and several PDE and variational calculus techniques have 
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been discovered to play also relevant roles to model several situations encountered in this context. 
We refer the reader to the nice book 0] for a quite comprehensive description of this field. One of 
the first effective contributions where variational methods and PDE were used specifically for the 
inpainting problem is due to Beltramio, Shapiro et al. 0. Several consequent papers came out and 
we refer to |11| and [3] as some of the most recent, interesting, and relevant for our purposes. 




Figure 2: The original template picture is illustrated on the left. Few color fragments are illustrated 
by darker disks on the right, where the gray level of the missing part is known. For the reader 
convenience, the gray level brightness has been exalted in order to distinguish it from the color part. 

All these important contributions are essentially addressed to solve the problem of guessing or 
learning some relatively small missing part of a signal/image by using the information of the relevant 
known part, and they are based on deterministic methods. The problem of learning from examples 
is a corner stone in information theory and artificial intelligence. Very recently a probability theory 
of what can be learned from a relatively large distribution of data according to a (unknown) proba- 
bility measure has been formalized in the beautiful contribution by Cucker and Smale 12 . Inspired 
by this work, subsequent papers, for example illustrate methods to construct estimators of the 
probabilistically best solutions of the learning problem. 

The situation just exemplified by the fresco problem is slightly different. 
We should reconstruct a large missing part, from a relatively small complete information and from a 
partial knowledge of the missing part itself. In fact, see for example Figure 1, what is known of our 
fresco is at least the gray level of the missing part. This additional information should be exploited 
to recover larger parts with respect to what was possible with the classical interpolation and PDE 
inpainting methods only. Moreover, one should expect that the additional information can increase 
the probability that good estimators have to reconstruct the signal/image, in the spirit of Cucker 
and Smale ideas. 

In this paper we want to show how interpolation and deterministic PDE inpainting methods 
can be separately adapted to work properly in this different situation, and how they can indeed 
interplay together to achieve nice results in concrete cases. Our aim here is not to give an exhaustive 
description of such technique, but mainly to introduce some new framework and models. Subsequent 
papers will detail theoretical and numerical properties of the methods here introduced. 

The paper is then organized as follows. Section 2 describes the general framework of the problem 
and its particular formulation in the case of the color reconstruction of digital images from few 
complete samples and a large information on the gray levels. Section 3 is devoted to illustrate a 
simple solution to the problem, based on interpolation techniques and nonuniform sampling methods. 
In Section 4 we discuss variational models and evolutionary PDE to solve the problem for ID signals 
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and for digital images. We conclude with the formulation of a system of coupled Hamilton- Jacobi 
equations for the solution of the color inpainting problem, its numerical implementation, and results. 
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Nota on color pictures 

This paper introduces methods to recover colors in digital images. Therefore the gray level printout 
of the manuscript does not allow to appreciate fully the quality of the techniques illustrated and 
the author recommends the interested reader to access to the electronic version with color pictures 
available online at [...] 

For the review process: Figure 15 and a corresponding animation are currently available online 
with colors at http://www.math.unipd.it/~mfornasi/abstractJMIV.html 

2 Scope of the problems 

The nonuniform sampling methods to recover a band-limited function u from an irregular set of 
sampling points S := {u(xj)}j^j described in |14l 1151 ITfil IB] are based on the assumption that the 
sampling points are dense enough with respect to the band-width of the function to be restored. This 
exactly in the same spirit as the more classical and well-known Whittaker-Shannon theorem. This 
is intuitively clear in the sense that if an image is rich of details and the missing part is quite large 
with respect to the known part, it should be impossible to recover correctly its original version. It 
is matter of the information that one has in fact at disposal. 




Figure 3: Voronoi decomposition of the domain where the nodes are the known complete color pixels. 
We are interested here to discuss and to treat a relatively different inverse problem. 
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Problem 1. Consider a signal u : SI — > K m , SI C K n , and a submanifold M. C K m of dimension 
d < m endowed with a C 1 nonlinear map £ : M m — > 7W. We call the couple (M, C) a C 1 -nonlinear 
projection. We assume to have a set of sampling nodes TV := C SI and a subset I C J 

such that C := {w(xi)}j e i, subset of R m , and I? := {£(u(xj))} j^jvr, subset of M, are known sets 
of values. The problem is then: How to reconstruct u from the sampling values CUV. 

Of course, the nonlinear distortion C is assumed as a datum of the problem, and in concrete cases 
one should have estimated it a priori. This is the typical situation where learning from examples 
|121 [5] is necessary, the nature of the real source of the distortion being unknown in general. In the 
following we will discuss briefly this estimation problem in the case of the color-gray conversion for 
images, exemplifying the process with the art fresco restoration (Figure 4). 

Examples 1. 1. Clearly if the function £| ran („) is injective, then it is invertible and one can 
construct a new set of sampling values S :— C U £u m / u \D. If S is dense enough and if u is a suitable 
band-limited function, then clearly it is possible to recover the function by the known nonuniform 
sampling methods. 

2. If G M. and C : W n — > M. has vanishing range ran(£| ran ( u )) = {0}, then the problem reduces 
again just to the nonuniform sampling problem where S := C. 

The interesting cases are when, for example, the function C is not injective on ran(tt), dim(A^) < 
m, and the density of the points {xi}i^x is not enough to ensure the perfect reconstruction of the 
function u. In relevant cases, even if the information contained in C is not sufficient to recover the 
function, maybe the additional information T> can be exploited "somehow" . 

This very general framework fits with many possible applications raising in concrete cases: For 
example, the recovery of a transmitted signal affected by a stationary (nonlinear) distortion, or, as 
in the case mentioned at the beginning of our paper, the restoration of colors of a fresco from its 
fragments and some gray picture taken prior to the damage. 

A digital image is modeled as a function u : SI C K 2 — > K+, so that, for each "point" x of the 
image, one associates the 3D vector u(x) = (r(x), g(x), 6(x)) G of the color divided between the 
different channels red, green, and blue. In particular a digitalization of the image u corresponds to 
its sampling on a regular lattice rl? . Let us again write u : M — > u(x) = (r(x), g(x), 6(x)), for 
x G W := SlflrZ 2 . Usually the gray level of an image can be described as a submanifold M C R 3 of 
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Figure 4: Estimation of the nonlinear curve L from a distribution of points with coordinates given 
by the linear combination ar + (3g + 76 of the (r, <?, b) color fragments (abscissa) and by the corre- 
sponding underlying gray level of the original photographs dated to 1920 (ordinate). The sensitivity 
parameters a, (5, 7 to the different frequencies of red, green, and blue have been chosen in order to 
minimize the total variance of the ordinates a 2 . Observe that, in this particular case, the curve L 
behaves nearly linearly in the interval [60, 170]. 

dimension 1 parameterized as M. := Ai T — {r(gl) : gl = £(r, g, b) := L(ar + /3g + jb), (r, g, b) G }, 
where a, f3, 7 > 0, a + /3 + 7 = 1, L : R — * K. is a non-negative increasing function, and t : M + — > M. 3 , 
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is a section such that C o t = idjj + . The function L is assumed smooth, nonlinear, and normally 
nonconvex and nonconcave. Therefore, in our terminology, (M,C) is a nonlinear projection, where 
C(r, g, b) :— L(ar + (3g + 76). For example, Figures 1 and 2 illustrate a typical situation where this 
model applies and Figure 4 describes the typical shape of an L function, which is here estimated by 
fitting a distribution of data from real color fragments. 

In fact, in this case, there is an area £l\D of the domain fl C M 2 of the image, where some 
fragments with colors are placed and complete information samples C are available, and an other 
area D where only a gray level information T> is known. 




Figure 5: The color of the pixel is extended in the direction of the corresponding Voronoi patch. 
The constraint check of the projection onto the sub-manifold up to a threshold allows to keep only 
the color information which is compatible with the known gray information. 

So our abstract Problem 1 can be reformulated as follows: When is possible to reconstruct 
the colors from the known information of the colors of the fragments C and from the gray level 
information T> of the missing part? There is not a unique answer to this question and different 
solutions are described depending on the methods chosen to solve the problem. Here we want to 
illustrate two techniques: The first based on interpolation methods and the second on PDE and 
variational calculus. Next we discuss how these two different strategies should be combined in order 
to achieve possibly better results. 




Figure 6: From the new deduced information, one applies again a Voronoi domain decomposition 
and iterates the process. 



3 Nonuniform sampling methods and interpolation techniques 

On one hand if the sampling set C is quite spread and dense, then clearly one might apply a 
nonuniform sampling method, maybe those described in [Tl 116) . or any other efficient interpolation 
method to recover the colors of the image. In this situation the information T> would not turn out 
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to be of any use. On the other hand, if the sampling set C is relatively small and concentrated, then 
the previous technique will fail or, at least, will give a very inaccurate result and therefore also the 
information T> should be taken into account and maybe exploited. 

Thus, we propose here a very simple scheme based on interpolation techniques that in the prac- 
tice gives acceptable results. Assume to have at our disposal the following procedures EXTEND, 
THRS, and ESTIM: 

EXTEND [F] — > G: given an incomplete color image F the procedure EXTEND computes a 
new color image G which coincides with F in the known color region of F and it extends the colors 
of F out of this region. 

THRS[F, F , C, e] — > G: given a color image F, a gray image F , a color-gray conversion func- 
tion C, and a threshold e > the procedure THRS computes a new image G which coincides with 
F only in those points x such that |£(F(x)) — Fo(x)| < e, and it is elsewhere. 

ESTIM[F, F ] — » {C,<7 2 }: given a color image F and a gray image F the procedure ESTIM 
estimates the color-gray conversion function C and the total variance of the ordinates a 2 with respect 
to the corresponding L function as described in Figure 4. 

At this moment we do not care of the accuracy the EXTEND procedure can achieve, nor its pos- 
sible concrete implementation. With these procedures at hand we define the following reconstruction 
algorithm: 

Algorithm 1. RESTORE[F, F ] -> G: 

G = 0; 

While G ^ F do 
G = F; 

{C,a 2 } = ESTIM[F,F ]; 

£ = r(a 2 ); 

F = EXTEND [F]; 

F = THRS[F,Fo,£,e]; 

od 

The procedure RESTORE[F, Fo] — * G computes an image G which is the color reconstruction 
of the image F from the samples C, being known also the gray level information T> of the missing 
part Fo. The threshold e > used at each iteration is a function of the variance a 2 , since we 
cannot expect to be more precise in our reconstruction than the intrinsic uncertainty of the model 
describing the color-gray conversion. 




Figure 7: Successive iterations of the interpolation process. The final image is illustrated on the 
right. Observe that neither the nose nor the eyes of the character could have been reconstructed. 
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In order to understand how the Algorithm 1 in fact works, let us illustrate here a simple imple- 
mentation of a procedure EXTEND. Clearly in more concrete applications one should use more 
sophisticated, accurate, and robust methods Given the color samples C of our image 

F one splits the domain by using a Voronoi decomposition, where the nodes for the decomposition 
are C, see Figure 3. One realizes a template procedure EXTEND just extending the color of any 
sample in C in the corresponding region of the Voronoi decomposition whose it is part, see Figure 
5. At this point we apply the procedure THRS just to keep the best possible extension which is 
compatible with the known gray level, up to the prescribed threshold. From this new extended set 
of color samples we iterate the process, constructing a new Voronoi decomposition, realizing the 
extension, and then applying again the thresholding and so on, see Figures 6 and 7. During the 
iterative scheme, a dynamical learning process of the model describing the color-gray conversion is 
realized. 

It is not difficult to see that Algorithm 1 converges after a finite number of iterations, and that it 
will stop having computed some possible extension G of the image F. This G might be considered 
the best possible extension of the samples C compatible with the known gray level samples T>. This 
method depends strongly on the way the procedure EXTEND is implemented and it exploits the 
information given by C and T> in an independent way. It is clear that one cannot reconstruct with 
this particular implementation any color which is not already included in the initial sample set C, 
and not connected geometrically with one of the initial sample color pixel by a sequence of Voronoi 
decompositions. Therefore, this method is highly local (i.e., no global properties of the image are in 
fact used) and it will fail whenever we will try to reconstruct colors which cannot depend somehow 
from those already known and given, and not connected geometrically with one of them. For this 
reason we want here to discuss also other possible inpainting techniques based on global properties 
of the signals. 



4 PDE and variational methods 

The modern techniques of processing digital signals are mainly and traditionally based on developing 
and applying Fourier and harmonic analysis concepts, e.g., wavelets |27| . time-frequency/Gabor 
analysis (171 1181 02] , sampling theory |51 [7] , and on stochastic and Bayesian modeling [28] . Only 
recently signal processing has become a very attractive field for applications of PDE and variational 
methods. We refer the reader to the nice introduction 0] for a presentation of this emerging field, 
for major details, and an extended literature. 

In this section we want first to show how Problem 1 can be modeled for ID signals by using 
suitable energy functionals whose minimum corresponds to possible solutions. In particular, we 
illustrate a numerical scheme which can be applied in the cases where the solution exhibits enough 
regularity. Next we propose a non-standard technique based on time-frequency/Gabor analysis to 
regularize problems involving non-smooth solutions. This method will require to move from a ID 
problem to an equivalent 2D problem. This will introduce us to the treatment of 2D digital signals 
(digital images) and to the formulation of a variational scheme for the solution of Problem 1 for 
digital images. 

Let us assume il C K, A, \l > 0, and C € C 1 (]R). We look for solutions u e W 1,p (Sl), for p > 1, 
of the following 

Problem 2. 

arginf \ F(u) = n [ \u(x) - u(x)\ 2 dx + A / \£(u(x)) - u(x)\ 2 dx + [ \Vu(x)\ p dx,u £ Vy 1,p (£l) \ 
{ Jn\D Jd Jn j 

(1) 

where u £ W l,p (Q) is the given observed signal, which is presumed to be correct on tt\D and 
distorted by the (nonlinear) function £ on D, see Figures 8,9. 
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Figure 8: As a nonlinear distortion we consider for example L(x) — 1.8(x + l) 2 (x — i) 2 . 




Figure 9: The continuous curve represents the original signal, while the dashed curve is the distorted 
signal by the nonlinear mapping. 

To minimize the functional F means essentially to find a function u £ W 1,p (£l) that coincides 
with u on Q\D, C(u) = u on D, and it minimizes some p-norm of its variation. In particular, 
for p — 2 we would look for solutions with small curvature. Let us fix p — 2 in the following for 
simplicity. By the direct method of the calculus of variations 4, Section 2.1] |13| it is well-known that 
Problem 2 has a solution u 6 W 1 ' 2 ^), and it solves the Euler-Lagrange equation 

Problem 3. 

= - Au + 2/i(u - u)l n \ D + 2\(£(u) - u) — (u)l D := E(£,u) (2) 

ax 

If C is convex then the solution of Problem 2 is unique, and Problem 3 is equivalent to Problem 
2. In the case C is not convex then we can anyway look for solutions of Problem 2 among those of 
Problem 3. 




Figure 10: The steepest descent method is applied to recover the original signal from the partial 
information on the missing part and the original information on the known part. In this nice example 
the reconstruction of the signal is rather good. 
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A way to compute a solution of Problem 3 is the steepest descent method: The solution u : 
f2 x R + — > K of the evolutionary equation 

^M = -£(A«(i,t)), u(.,0)=«o. (3) 

tends to make vanish £(C,u(x,t)) for t — > +00. Therefore one can look for solutions of Problem 2 
given by u(x) :— lim^oo u(x, t). 

REMARK: This method is as more effective as uq, the initial guess, is closer to the real solution 
u. As we will discuss later, this is the crucial and important motivation to combine interpolation 
methods, for achieving the best first guessing of the solution, and then variational techniques to 
complete the reconstruction. In fact the result shown in Figure 10 would have been impossible by 
using an interpolation (or irregular sampling) technique only, since the size of the gap (missing part) 
to be reconstructed would have been too large with respect to the band-width of the signal, and 
the use of the variational method might fail if the first guess uq is not close enough to the solution. 
For example, in Figure 11 we show the evolution of u(x,t) in reconstructing the missing part of 
the signal in Figure 9, from a first linear interpolation guessing. Therefore, we expect that in con- 
crete and practical problems only an interlacing of these two different tools can give efficient results. 

The steepest descent method illustrated in Q can be easily discretized and it can work on 
samples of the signal by means, e.g., of an explicit Euler scheme. Let us denote A 2 (v)(i) := 
v(i — 1) — 2v(i) + v(i + 1), for any vector v and for i = 1, N — 2 and define the procedure 

£[£, v, i] A 2 (v)(i) - 2 M (v(i) - u(i))l n \ D - 2A(£(v(i)) - u(i))^(v(i))l D , i = 1, N - 2. 

Algorithm 2. STEEP_DESC[u , e, C, At] : 

v = v := u ; 

While max{£[£,v,i]|,i= l,...,A-2} > e do 
For i=l,...,N-2 do 

v(i)=v(i) + Ai (£[Av,i]); 

od 
v = v; 

od 

U = V. 

Therefore a solution of our original Problem 1 can be modeled as the result of the procedure 
STEEP_DESC[uo, £, C, At] which, given an observed set of N samples Uo = u := C U T>, where C 
are presumed correct and T> are presumed distorted by the function C, computes approximate evo- 
lutions of Uo by discrete time steps At. Once the Algorithm 2 has corrected the distorted samples, 
one can use any interpolation result to reconstruct the solution u of the Problem ^ 

A similar approach can be realized also in the case the solution is not smooth, for example 
u G BV(£l). In general the discretization method to compute a minimizing solution u S BV(fl) of 
the Problem^is more sophisticated and usually it is necessary to construct an associated regularized 
problem which has minimum approximating u. We refer, for example, to 01 Section 3.2] for details 
on techniques used in this situation. 

In the next section we want to present a non-standard regularization method based on time- 
frequency and harmonic analysis tools. 
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Figure 11: Iterations of the discretized steepest descent method 



4.1 A non-standard regularization method: time- frequency analysis 

For a function / e L 2 (R n ) and 0^.g£ C£°(R"), \\g\\ 2 = 1, the transformation 



V g (f)(x,u):= f f(t) e - 2 ^ t g(t-x)dt, 



(4) 



is called the short time Fourier transform (or windowed Fourier transform or simply Gabor trans- 
form) of the function /. Since g is a smooth function, even if / is not smooth, V g (f) is smooth 
indeed. Moreover, the map V g : L 2 (M. n ) — > L 2 (R 2 ™) is a unitary isomorphism into its range and 
lcft-invertible by its adjoint 



V* (F) = [ F(x, u)e 27!lut g{t - x)dxdu>. 

JM 2 " 



(5) 
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Figure 12: Iterations of the discretization of the steepest descent method in the time- frequency plane 



Therefore, one has the following reproducing formula 

f(t) = V g * o V g (f)(t) = f V g (f)(x,u>)e M g(t - x)dxdw. 



(6) 



Moreover II := V g oV* is the orthogonal projection from L (R ) onto ran(V g ). The Gabor transform 
gives simultaneous information on the time/space- frequency content of a given signal. In particu- 
lar, it tells us which are the "instantaneous" frequencies appearing in a ID signal at a given time, 
somehow as the score is describing a music, telling which are the notes to be played at a given time. 

Because of the invertibility of V g , to analyze or synthetize a signal in the time, or to do that in 
the time-frequency plane are equivalent operations. For example, in Figure 12 we have visualized the 
Gabor transform of the iterations in Figure 11. This suggests that Problem [21 can be reformulated 
in the time-frequency plane. In fact, as already mentioned, it is known that u £ L 2 (fl) implies 
V g (u) € L 2 (£l + supp(g) x R) fl C°° . The treatment of the problem as defined in the time- frequency 
plane is non-standard, especially because the domain itself is unbounded. On the other hand, the 
range of V g is always constituted by smooth functions and this of course is a way to regularize the 
problem. 



Problem 4. 



arginf{F(C/) = (i [ 
Jn 



n g \D g xi 



\U(x,(jj) — V g (u)(x, ui)\ 2 dxduj+\ / \Ctf(U)(x,uj) — Vg(u)(x,oj)\ 2 dxduj 



D> 



\VU\ 2 dxduj,U € W 1,2 (il g 



(7) 
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where n g = fi + supp(g), D g = D + supp(.g), Ctf{U) := V g (C(V*(U))), and u e L 2 (fl) is the given 
observed signal, which is presumed to be correct on Q\D and distorted by the (nonlinear) function 
C on D. 

As usual, one can correctly compute the associated Euler-Lagrange equations and to formulate 
the corresponding steepest descent method, in the same way as we have done in the ID case. Since 
similar definitions and properties of the Gabor transform can be formulated for discrete signals one 
can easily discrctize the steepest descent method and define a suitable numerical scheme. We refer 
the reader to ^| ^] for major details on numerical Gabor analysis in the treatment of digital 
signals. The analysis of Problem 01 and its connections with corresponding properties in the time 
domain will be investigated in successive contributions. 

This section has been useful to us also to move our attention from the ID situation back to a 2D 
problem, which was the inspiration of this paper. In the next section we want to generalize what 
we have discussed for univariate signals on the real line to multivariate signals on 2D domains. In 
particular, it is well known that the Laplacian appearing in the Euler-Lagrange equations associated 
to Problem 0] has very strong isotropic smoothing properties and does not preserve edges. While 
this could not be a relevant problem for 2D function in the range of V g , since they are smooth 
and affected by the Heisenberg uncertainty principle which makes them intrinsically "blurred" , it is 
not suitable for the reconstruction of 2D functions representing natural images, which are usually 
characterized by the presence of edges, curves, and maybe fractal structures. 



4.2 Variational models for image inpainting 

There are several and different variational methods for solving the so called image inpainting problem, 
i.e., the reconstruction of a small missing part of an image by using the information of the remaining 
relevant known part. Each of them offers some nice properties and effects together with some 
drawbacks. In this section we want to propose a modification of what we consider the most simple 
and well known approach, in order to have a reasonable numerical solver for our original problem 
of the reconstruction of colors. First, let us describe in the following three of the models discussed 
recently in the literature: 

1. The Beltramio-Shapiro-Caselles-Ballester model. Let L(u) be a smoothness measure of the 
image u, for example, a second order differential operator given by 

L(u) := /(Vu,V<g> Vu). (8) 



Typical choice for L is the Laplacian L(u) = Au = trace(V <8> Vu). The reconstruction model 
proposed by Beltramio, Shapiro et al. [H] is based on the evolutionary equation 

du{{X dt y) ^ = VVfoy),*) • VL(u(x,y),t)), (9) 

where V^u = (- Mg-lM , = |Vu|T, point to the tangent T. This means that, 

at the equilibrium, one has dI g^ = 0. Thus, in terms of boundary smoothness data, the 
inpainting process evolves transporting the information along the extended isophotes (i.e., level 
curves)into the region to be reconstructed. This method has a nice numerical implementation 
in [Hj which allowed the authors to show several examples where it works nicely. However, 
since there is not a clear information connecting different isophotes, the evolution might create 
nonexistent T-junctions inside the inpainting region and smoothing operations are necessary 
to connect properly such discontinuities. Moreover, a rigorous theoretical description of such 
equation is still a matter of investigation. 
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The Chan and Shen non-linear diffusion model. This model is much more in the spirit of the 
energy methods presented in the ID model and it offers several advantages. In particular, the 
mathematics related to such model is much more developed and numerical implementations 
have been recently well formulated. It is based on the following minimizing problem 



arg inf 




Q\D 



u{x) -u(x)\ 2 dx + / ct>(\Vu(x)\)dx,u£W 1 ' p (n) 



(10) 



where u S BV(Q) is the given observed signal, which is presumed to be correct on Q\D and 
D is the region to be reconstructed. 

It is known that a solution of the problem (see for example Proposition 3.3.4, Section 3.3.1]) 
that a minimizing solution solves in particular the following Euler-Lagrange equation 

>'(|Vu| 







|Vu| 



Vu + 2{u{x) - u(x))l QV 



(11) 



endowed with suitable boundary conditions, where Vu is the approximate derivative of u. 

On one hand, such an evolutionary equation and the corresponding steepest descent method 
are normally used for the numerical solution of the energy problem as reasonable practical 
methods. On the other hand, this procedure is not completely mathematically correct and 
more sophisticated techniques should be implemented and we refer the reader to @] Section 
3.2.3-4] for more details. Such techniques are based again on regularization methods. It has 
been also proposed by Chan and Shen to modify equation Ijllf) in the following way 

(|Vu|), 



= - Vw V- 



IVuI 



Vu +2(u(aO-u(aO)l n \B. 



(12) 



The steepest descent method in the inpainting region D, for the trivial choice <j)(t) = t, becomes 



du 

at 



— = iVitlV 



Vu 



(13) 



which describes the well-known (morphological invariant) mean curvature motion, see, e.g., 
1101 II 1 j . Efficient numerical implementations of such an evolutionary equation have been 
recently formulated for example in j2Ej. Such model tends to approximate the isophotes in the 
inpainting region by straight lines and, in the following, we will show that one can modify it 
to avoid this problem. 

The Ambrosio and Masnou model. A method to extend in a smooth way the isophotes from 
T-junctions at the border of the inpainting region has been recently proposed and studied by 
Ambrosio and Masnou .JJJ based on the following minimization problem 



arg inf If(u) = J_\Vu\(l 



Vu 



dx,ue W^ p (9) 



(14) 



where D is a slightly larger domain containing the inpainting region D. In fact, in their paper 
[3] they discuss the corresponding relaxed functional F with respect to the L 1 convergence and 
the existence of BV(Cl) minimizers of F. In the case n = 2 and p > 1 Chan and Shen derived 
the (fourth order) Euler-Lagrange equation corresponding to such energy problem |llj . 

From what we have just illustrated above, there are several nice intuitions and models based on vari- 
ational methods in order to solve the inpainting problem, where deep mathematics and numerical 
problems are involved and still rather unexplored. In our understanding, presently the mathemat- 
ically most well supported and the most promising in terms of efficient numerical schemes and 
realizations is the Chan-Shen approach. We want now to borrow their approach in order to discuss 
and to model our color reconstruction problem. We assume then that u is a RGB image. 
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Problem 5. 



arginf lF(u) = fi / \u{x) - u(x)\ 2 dx + A / \£{u{x)) ~ u{x)\ 2 dx + / 0(|V«(x)|)da;,« G W^fSl,^) 

(15) 

where u is the given observed image, which is presumed to be colored on £l\D and just gray level 
on D where the transformation from RGB to gray level is given by the (nonlinear) function C. 

Since we are dealing with natural images one should relax the functional and look for BV min- 
imizers. This strategy can be handled in a similar way as it is done, for example, in Theorems 
3.2.1-2] and it will be discussed elsewhere. However it is true that the BV solution of the relaxed 
formulation of Problem [S] is anyway solving the following Euler-Lagrange equations 



k = 1,2,3, 



= -V- ( ^f|^"| ^ v "fc) + 2 ^fc ~ u k )ln\D + 2\{£(u) - u)^-(u)l D :=£ k (£,u), 

. (16) 

where u\ :— r, ui := g, U3 :— b are the RGB components of the image u, and Vw^ is the approximate 
derivative. One can consequently formulate the three evolutionary steepest descent equations by 

dUk{i ^ v),t) = -£ k (jr.,u((*,y)M u(;0) = u ,k = 1,2,3, (17) 

to make vanish £k(£-,u((x,y),t)), k — 1,2,3, for t — > +00. Therefore one can look for solutions of 
Problem [S] given by u{x,y) := lim^oo u{{x, y), t). Observe that, in particular, at the equilibrium 
and on the inpainting region D one has 

= -V • (^pV« fc ) + 2A(£(«) - u)^-(u). 

This means that the extension of the isophotes will not be in general just straight lines, as we have 
already seen also in the ID model (Figure 10) where the reconstructed part is not just a linear 
interpolation as one could expect using a mean curvature motion evolutionary equation. 

It is important to observe that (|16fl and (|17fl are systems of coupled second order Hamilton- Jacobi 
equations and the analysis of the solutions constitutes itself an open problem of independent inter- 
est that will be discussed elsewhere. Some results on existence and uniqueness of viscosity solutions 
for systems of steady nonlinear second order PDE have been proposed by Koike |24l 1251 . Un- 
fortunately those results apply under quasi-monotonicity conditions that here could not in general 
be verified. However, one can easily construct solutions of i|16[) and l|17|) in particular cases as follows. 

Let us assume that a gray level image v can be identified with an RGB image r o v, by means 
of a section r : K + — > ]R^_ such that C o r = id. For example, if £(r, g, b) = ar + fig + 76, with 
a + (3 + 7 = 1 then one can choose r(gl) = (gl, gl, gl). Suppose that u is an image with piecewise 
linear/straight isophotes, colored on Q\D and gray on D. Therefore we can always assume that 
u = (iti, 1*2, W3) = U|o\d + t { u \d) is a vector valued function. For the choice <p(t) — t, it is not 
difficult to show that such a function u is almost everywhere a stationary solution of l|17|) . i-e., for 
the choice uq = u := u there will be no evolution, and u is solution of Ijl6|l . 

4.3 The numerical implementation 

The steepest descent method illustrated in Ijl7|l can be easily discretized and it can work on samples 
of the image by means, e.g., of an explicit Euler scheme. We assume here 4>{t) = t and we modify 
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Figure 13: Successive iterations of the nonlinear diffusion process. Starting from Figure 2 the color is 
diffused extending the isophotes in the inpainting region, following the tracks of the underlying gray 
level. Emblematic is the reconstruction of the color on the nose of the character, where the diffusion 
progresses along it, but not through its boundary, the isophotes having there almost vanishing 
curvature. The last image on the bottom-right is the result of 300 iterations of the outer loop of 
Algorithm 3. Here we have chosen At = 0.1, A = /z = 10, and £(r,g,b) = j(r + b + g). The initial 
information content does not appear sufficient for a correct color reconstruction and only a limited 
part of the image can be partially restored. 



again the velocity of the evolution as in (|12|) by multiplication of the first term with the gradient. 
This produces a motion of mean curvature of the isophotes of the images. 

For a color digital image v = (vi, V2, V3) and for i = 1, N — 2, j — 1, M — 2, let us denote 
/ \ / • • \ 2(v fc (m,n) -v k (i,j)) na s 

( * K ' J) ■=„,„£<.„> IVv.^lMVv^,, < 18 > 

where Af(i,j) denoted the 4-neighborhood of the pixel position and Vu(i, j)\ — (S^.u(i, j)) 2 + (S^u(i, j)) 2 

S 1 ^ , Sy being the central differences of order 1 in the direction x and y respectively. Then define the 
procedure 

£ k [£,v,i,j] ->■ |Vv k (i,j)|«;(v k )(i,j) - 2/i(v k (i,j) - u k (i,j))l \ D - 2A(£(v(i,j)) - u(i,j))- — (v(i,j))l D , 

dx k 

i = l,...,N-2, j = l,...,M-2, k= 1,2,3. 



1G 




Algorithm 3. STEEP_DESC_2D[u , s, C, At] : 

v := v := u ; 

While max{|£[£,v,i,j]|,i= l,...,JV-2,j = l,...,M-2} > e do 
For i=l,...,N-2 do 
For j=l,...,M-2 do 
For k=l,...,3 do 

v*(i, j) = v fc (i,j) + Ai (S k [£,v,i,j]); 

od 

v = v; 

od 

od 

od 

U = V. 




Figure 14: Successive iterations of the nonlinear diffusion process. Starting from a few colored 
samples more the diffusion process is able almost to complete the reconstruction of the color of all 
the face of the character after 300 iterations. Again we have chosen here At = 0.1, A = \i = 10, and 
C(r,g,b) = \{r + b + g). 

Therefore a solution of our original Problem 1 can be modeled as the result of the procedure 
STEEP_DESC_2D[u , e, C, At] which, given an observed set of N x M samples u =Q:=CUD, 
where C are presumed colored and T> are presumed gray, computes approximate evolutions of Uo by 
discrete time steps At. Of course the explicit Euler scheme is not the most efficient, and implicit 
variants can be easily derived, see, e.g., |26|. 
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Let us show now some applications ol Algorithm 3 for the color inpainting problem in Figure 13 
and Figure 14. 

4.4 Conclusion: The interpolation-inpainting method 

Our reconstruction method is based on an inpainting model which exploits additional information 
given by a known nonlinear projection of the solution onto a manifold of lower dimension. The 
numerical solution can be implemented by a suitable discretization of evolutionary equations which 
are assumed to start from an initial guess of the solution. We have shown that in concrete problems 
where the inpainting region is large with respect to the known complete data, only a combination of 
suitable interpolation and variational methods can work properly. Since the interpolation methods 
can be computationally expensive (for example the construction of a Voronoi decomposition is indeed 
rather expensive) one can combine a few initial iterations of Algorithm 1, in order to extend as much 
as possible the knowledge about the color-gray conversion and a possibly good initial guess of the 
solution, with successive fast iterations of Algorithm 3 to complete the reconstruction, see Figure 
15. 




Figure 15: Successive iterations of the interpolation-inpainting process. The first three pictures on 
the top are the result of iterations of Algorithm 1, and the successive are generated by Algorithm 3. 
Again we have chosen here At = 0.1, A = \i = 10, and C(r, g,b) = |(r + b + g). 
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